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Abstract 

We present a measurement of the underground decoherence function using multi-muon 
events observed in the MACRO detector at Gran Sasso at an average depth of 3800 hg/cm 2 . 
Muon pair separations up to 70 m have been measured, corresponding to parent mesons 
with Pj_ < l-=-2 GeV/c. Improved selection criteria are used to reduce detector effects 
mainly in the low distance separation region of muon pairs. Special care is given to a new 
unfolding procedure designed to minimize systematic errors in the numerical algorithm. The 
accuracy of the measurement is such that the possible contribution of rare processes, such 
as ^ + N — > ^ + N + /i + + can be experimentally studied. 

The measured decoherence function is compared with the predictions of the hadronic 
interaction model of the HEMAS Monte Carlo code. Good agreement is obtained. We 
interpret this agreement to indicate that no anomalous P± components in soft hadron- 
Nucleus and Nucleus-Nucleus collisions are required by the MACRO experimental data. 
Preliminary comparisons with other Monte Carlo codes point out that the uncertainties 
associated with the hadronic interaction model may be as large as 20%, depending on the 
energy. MACRO data can be used as a benchmark for future work on the discrimination of 
shower models in the primary energy region around and below the knee of the spectrum. 

I. INTRODUCTION 

The knowledge of hadronic interaction processes plays a fundamental role in studies of 
cosmic rays in the VHE-UHE range (10 12 eV < E < 10 17 eV). In particular, the interpre- 
tation of indirect measurements intended to determine the features of primary cosmic rays, 
such as spectra and composition, depends on the choice of the hadronic interaction model 
adopted in the description of the atmospheric shower development. For instance, muons 
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observed by deep underground experiments are the decay products of mesons originating 
mostly in kinematic regions (high rapidity and high yfs) not completely covered by existing 
collider data. The problem is particularly important for nucleus-nucleus interactions for 
which available data extend only to a few hundreds of GeV in the laboratory frame. It is 
therefore crucial to find physical observables which are primarily sensitive to the assumed 
interaction model rather than to the energy spectra and chemical composition of primary 
cosmic rays. 

The shape of the muon lateral distribution is well-suited for this purpose. In particular 
it allows the study of the transverse structure of hadronic interactions, which is one of the 
most relevant sources of uncertainties in the models 0. In fact, different aspects of the 
interactions contribute to the lateral distribution. We can qualitatively understand this by 
simple arguments, valid in a first order approximation. Let us consider a single interaction 
of a primary nucleon of total energy Eo, producing mesons of energy E n,K with transverse 
momentum P±, at a slant height H pro d, which eventually decay into muons. Calling r the 
separation of a high energy muon (i.e. moving along a straight line) from the shower axis, 
we have: 



In this simplified description we are neglecting the transverse momentum in the parent decay. 
The previous expression can be written in a more instructive way, considering that at high 
energy, apart from terms of the order of (mx/En) 2 , the longitudinal cm. variable Xp is 
approximately equal to the laboratory energy fraction: 



The assumption of an exponential atmosphere has been used in the last expression. It 
can be seen how the transverse and longitudinal components of the interaction, as well as 
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the inclusive and total cross sections, convolve together (with different weights) to yield 
the lateral separation. The role of P± remains a dominant one in determining the relative 
separation of the muon component by introducing a loss of collinearity ( "decoherence" ) with 
respect to the direction of the shower axis. 

A qualitative extension to the case of nuclear projectiles can be made within the frame- 
work of the superposition model, where each nucleon of the projectile of mass number A 
is assumed to interact independently with energy Eq/A. Further refinements are needed to 
account for modifications in the P± and xp distributions deriving from the nuclear struc- 
ture of projectile and target, as will be discussed later. A reliable evaluation of the lateral 
distribution function can be obtained only by Monte Carlo methods. 

Deep underground experiments are capable of selecting atmospheric muons in the TeV 
range produced in the initial stages of the extensive air shower (EAS) development. They 
can perform a measurement of muon separation which is highly correlated to the lateral dis- 
tribution. Since the shower axis position is not usually known, the distribution of muon pair 
separation in multimuon events is studied. Muons associated with the same events, coming 
in general from different parent and shower generations, are grouped together. Further- 
more, a wide range of primary energy is integrated in the same distribution. It is generally 
assumed, and supported by many simulations, that the shape of this distribution is only 
slightly affected by the mass composition of primaries 0, thus preserving the sensitivity to 
the interaction features. As an example, in Fig. [l] we show the dependence of the average 
pair separation, as detected at the depth of the underground Gran Sasso laboratory, with 
respect to the (P±) of the parent mesons and to their production slant height in the at- 
mosphere. These have been calculated by means of the HEMAS Monte Carlo code |J for 
a mixed primary composition [|l(J. This code employs an interaction model based on the 
results of the experiments at hadron colliders. 

The decoherence function as measured in an underground experiment is also affected by 
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multiple scattering in the rock and, to some extent, geomagnetic deflection. 

For a detector with geometrical acceptance A(8,(f)), for zenith and azimuthal angles 9 
and 0, respectively, we define the decoherence function as the distribution of the distance 
between muon pairs in a bundle: 

d N_ 1 t 1 



dD VtT J A(9,<j>) dDdVL 

where N(D,9,4>) is the number of muon pairs with a separation D in the direction (9,<p), 
Q is the total solid angle covered by the apparatus and T is the total exposure time of 
the experiment. A muon bundle event of multiplicity iV M will contribute with a number of 
independent pairs N = N^N^ — l)/2. 

In principle, a decoherence study can be performed without a single large area detector, 
and in early attempts the muon lateral separation was studied via coincidences between two 
separate movable detectors ||. The advantage offered by a single large area detector is the 
ability to study the features inherent in the same multi-muon event, such as higher order 
moments of the decoherence distribution ||. 

The large area MACRO detector || has horizontal surface area of ~ 1000 m 2 at an 
average depth of 3800 hg/cm 2 of standard rock (E^ > 1.3 TeV) and is naturally suited 
for this kind of measurement. An analysis of the muon decoherence has already been per- 
formed [p|,|5|,P|,|TU[ . The bulk of multiple muon events in MACRO corresponds to a selection 
of primary energies between a few tens to a few thousands of TeV/nucleon. Hadronic in- 
teractions and shower development in the atmosphere were simulated with the previously 
noted HEMAS code. In particular, a weak dependence on primary mass composition was 



confirmed for two extreme cases: the "heavy" and "light" composition models [11]. The 
MACRO analysis was designed to unfold the true muon decoherence function from the 
measured one by properly considering the geometrical containment and track resolution ef- 
ficiencies. This procedure permits a direct comparison between measurements performed by 



different detectors at the same depth, and, more importantly, whenever new Monte Carlo 
simulations are available, allows a fast comparison between predictions and data without 
the need to reproduce all the details of detector response. 

The first attempt, obtained while the detector was still under construction, and there- 
fore with a limited size, was presented in . The same analysis, with a larger sample based 
on the full lower detector, was extended in M. With respect to the HEMAS Monte Carlo 
expectations, these results indicated a possible excess in the observed distribution at large 
separations. In Ref. [nj we presented the decoherence distribution without the unfolding 



procedure; the claimed excesses were not confirmed. In order to reach more definitive con- 
clusions, a more careful analysis of the systematics associated with the unfolding procedure 
was considered necessary. A detailed discussion of this item will be addressed in Section 4. 

A more careful discussion of the Monte Carlo simulation is also necessary. The bulk 
of the muon bundles collected by MACRO are low multiplicity events, coming from parent 
mesons in the far forward region of UHE interactions, not easily accessible with collider 
experiments. This requires an extrapolation to the highest energies and rapidity regions, 
introducing possible systematic uncertainties. For instance, some doubts have been raised 
ll concerning the treatment of meson P± in HEMAS. In the HEMAS hadronic interaction 
code, secondary particle P± depends upon three different contributions: 

• (P_l) increases with energy, as required by collider data in the central region; 

• (P±) increases in p-Nucleus and Nucleus-Nucleus interactions, relative to that for pp 
collisions, according to the "Cronin effect" |jl2| ; 



(P±) varies with xf, according to the so called "seagull effect" [p~3[ . 



The sum of these effects yields some doubt about a possible overestimate of P± for energetic 
secondary particles, an hypothesis recently restated in WM. It is therefore crucial to perform 
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a high precision test of the transverse structure of this model, since it affects the calculation of 
containment probability for multiple muon events and, consequently, the analysis of primary 
cosmic ray composition 

In this paper, a new analysis of the unfolded decoherence function is presented, performed 
with improved methods up to 70 m. The present work enlarges and completes the data 
analysis presented in [§,[!(]]. Preliminary results of this unfolding procedure showed an 



improved agreement between experimental data and Monte Carlo predictions. 

Particular attention is paid to the small-separation (D < 1 m) region of the decoherence 
curve, in which processes such as muon-induced hadron production can produce a back- 
ground to the high energy muon analysis. At the energies involved in the present analysis 
{E^ > 1 TeV), moreover, muon-induced muon pair production in the rock overburden could 
yield an excess of events with small separation, as suggested in ||16|1 . This process is usu- 



ally neglected in Monte Carlo models commonly adopted for high energy muon transport 



Section 2 is devoted to the description of the detector and of data analysis, with a focus 
on new event selection criteria. In Section 3 the features of the Monte Carlo simulation are 
presented together with the comparison between experimental and simulated data in the 
MACRO detector, while Section 4 is dedicated to the unfolding procedure. A comparative 
discussion of the features of different hadronic interaction models is summarized in Section 5. 
In Section 6, the problem existing in the first bins of the decoherence distribution is presented 
in detail, testing new hypotheses on its origin. Conclusions follow in Section 7. 

II. DETECTOR DESCRIPTION AND DATA ANALYSIS 

The MACRO detector ||, located in hall B of the Gran Sasso Laboratory, is a large area 
detector equipped with streamer tube chambers, liquid scintillation counters and nuclear 
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track detectors arranged in a modular structure of six "supermodules" . Each of these is 12 
mxl2 mx9 m in size and consists of a 4.8 m high lower level and a 4.2 m upper "attico". 
In this paper only data from the lower level of the apparatus are included; therefore only 
the lower detector will be described further. 

Tracking is performed by means of limited streamer tubes, which are distributed in ten 
horizontal planes separated by ~ 60 g cm~ 2 of CaC03 (limestone rock) absorber, and in six 
planes along each vertical wall. The streamer tubes have a square cross section of 3x3 cm 2 , 
and are 12 m long. From each plane two coordinates are provided, the wire (perpendicular 
to the long detector dimension) and strip views. The latter employs 3 cm wide aluminum 
strips at 26.5° to the wire view. The average efficiencies of the streamer tube and strip 
systems were 94.9% and 88.2% respectively, in the period of this analysis. 

The spatial resolution achieved with this configuration depends on the granularity of 
the projective views. The average width of a cluster, defined as a group of contiguous 
muon "hits," is 4.5 cm and 8.96 cm for the wire and strip views, respectively. Muon track 
recognition is performed by an algorithm which requires a minimum number of aligned 
clusters (usually 4) through which a straight line is fit. The differences between the cluster 
centers and the fit determine a spatial resolution of <Jw=\.\ cm for the wire view and 
(75=1.6 cm for the strip view. These resolutions correspond to an intrinsic angular resolution 
of 0.2° for tracks crossing ten horizontal planes. 

In reconstructing the best bundle configuration, the tracking package flags track pairs 
as parallel, overlapping, or independent and not parallel. This is achieved in two steps, in 
each projective view: 

• two tracks are defined as parallel if their slopes coincide within 2 a or if their angular 
separation is less than 3° (6° if the tracks contain clusters whose widths exceed 30 cm). 
Otherwise, the track pair is flagged as independent and not parallel if its distance 
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separation is larger than 100 cm. 

• tracks at short relative distance are labelled as overlapping if their intercepts with the 
detector bottom level coincide within 3.2 a (2 a if their angular separation is < 1.5°). 

The routine chooses the most likely bundle as the set having the largest number of parallel 
tracks and the largest number of points per track. Subsequently, tracks flagged as not 
parallel are considered in order to include fake muon tracks originated primarily by hadrons 
or 5-rays in the surrounding rock or inside the detector. A two-track separation of the order 
of 5 cm is achieved on each projective view. However, this capability can be substantially 
worsened in case of very large, but rare, catastrophic energy losses of muons in the detector. 

Only tracks with a unique association in the two views can be reconstructed in three 
dimensional space. At this level, pattern recognition is used to require a complete matching 
between tracks belonging to different projective views. This is automatically achieved when 
two tracks pass through separate detector modules. When they are in the same module, 
matching of hit wires and strips on the same detector plane is accomplished by taking 
advantage of the stereo angle of the strips with respect to the wires. In some cases the 
track pattern correspondence between the two views is also used. The possibility to analyse 
muon decoherence in three dimensional space is important to have an unbiased decoherence 
distribution. However, the unambiguous association of muon tracks from the two projective 
views cannot be accomplished for high multiplicity events because, in events characterized 
by a high muon density, the tracking algorithm is not able to resolve the real muon pattern 
without ambiguities, especially when tracks are superimposed. In Ref. f9,[TC] we presented 



the muon decoherence function in the wire view alone, which allowed the extension of the 
analysis to higher multiplicities. 

We have analyzed about 3.4 • 10 5 events, corresponding to a 7732 hr live time for the 
lower part of the apparatus. These events were submitted to the following selection criteria: 
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1. Zenith angle smaller than 60°. This choice is dictated by our limited knowledge of 
the Gran Sasso topographical map for high zenith angles. Moreover, we cannot disre- 
gard the atmosphere's curvature for larger zenith angles, which at present our current 
simulation models do not include. 

2. Fewer than 45 streamer tube hits out of track. This selection is designed to eliminate 
possible misleading track reconstruction in events produced by noise in the streamer 
tube system and/or electromagnetic interactions in or near the apparatus. 

3. Track pairs must survive the parallelism cut. This rejects hadrons from photonuclear 
interactions close to the detector, as well as tracks reconstructed from electromagnetic 
interactions which survived the previous cut. 

The last cut is not completely efficient in rejecting muon tracks originating from local 
particle production because the angle between these tracks may fall within the limits imposed 
by the parallelism cut. These limits cannot be further reduced since the average angular 
divergence due to multiple muon scattering in the rock overburden is about 1° at the MACRO 
depth. This is a crucial point, since these events could contaminate the decoherence curve 
in the low separation region and are not present in the simulated data because of the 
excessive CPU time required to follow individual secondary particles. A similar effect could 
be produced by single muon tracks with large clusters, which may be reconstructed as a 
di-muon event by the tracking algorithm. 

In order to reduce these effects, a further selection was applied. We computed, for each 
muon track in the wire view, the ratio R between the number of streamer tube planes hit by 
the muon to the number of planes expected to be hit considering the track direction. Only 
tracks with R > 0.75 were accepted. The application of this cut (hereafter cut C4) in the 
wire view alone is a good compromise between the rejection capability of the algorithm and 
the loss of events due to the unavoidable inefficiency of the streamer tube system. We found 
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that in the wire projective view the probability to reject a muon track due to contiguous, 
inefficient planes is 2.0%. 

To show the effects produced by cut C4, we present in Fig. |] the fractional differences 
between the experimental decoherence curve before and after its application. As expected, 
the new cut affects only the first bins of the distribution. 

To test the ability of cut C4 to reject hadronic tracks, we used FLUKA |17j to simulate 
3028 hr of live time in which muons were accompanied by hadronic products of photonuclear 
interactions in the 10 m of rock surrounding the detector. We found that the parallelism cut 
alone provides a rejection efficiency of about 54.6% of the pair sample, while the addition of 
cut C4 enhances the rejection to 95.9%. The effect of hadron contamination, furthermore, 
is very small, contributing less than 1% in the overall muon pair sample. This estimate, 
together with the plot of Fig. |2| and the results of a visual scan, suggest that the main track 
sample rejected by cut C4 is made of large cluster tracks. After the overall application of 
these cuts, the number of surviving unambiguously associated muon pair tracks is 355,795. In 
Fig. ^| the percentage of the reconstructed events as a function of muon multiplicity is shown 
(open circles). In the same figure, the percentage of the unambiguously associated muon 
pairs as a function of the multiplicity is also reported (black circles). Due to detector effects, 
the number of associated pairs N' ir in an event of multiplicity is generally smaller than 
the maximum number of independent pairs N pa i r = N^N^ — l)/2. This reduction becomes 
greater for high multiplicities, for obvious reasons of track shadowing. In any case, we 
still find that the weight of high multiplicity events remains dominant in the decoherence 
distribution. In order to reduce this effect and to reduce the possible dependence on primary 
composition, we have assigned a weight 1/N' ir to each entry of the separation distribution. 
This prescription, followed also for simulated data, has been already applied in most of the 
previous analyses performed by MACRO. Moreover, we emphasize that the focus of this 
analysis is centered on the shape of the distribution; the absolute rate of pairs as a function 
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of their separation is neglected. 



III. MONTE CARLO SIMULATION 

The Monte Carlo chain of programs used in the simulation consists of an event generator, 
capable of following the development of the hadronic shower in the atmosphere and the muon 
transport code in the rock overburden, and a detector simulation package. We have used, as 
in all previous relevant analyses of muon events in MACRO PJIUfl, the HEMAS code |J as 
an interaction model and shower simulator. Nuclear projectiles are handled by interfacing 
HEMAS with the "semi-superposition" model of the NUCLIB library [[HJ. The final relevant 
piece of simulation is the three-dimensional description of muon transport in the rock. A 



comparison of the performance of different transport codes, reported in Ref . |2l| , showed that 
the original package contained in the HEMAS code was too simplified, leading, for instance, 
to an underestimated muon survival probability at TeV energies. In order to verify possible 
systematics affecting the decoherence distribution, we repeated the Monte Carlo production, 
interfacing the more refined PROPMU code |nj to HEMAS. We have verified that, at least 
to first approximation, no changes in the shape of the decoherence function are noticeable 
between the two different simulation samples. For this reason, the sum of the two different 
Monte Carlo productions will be used in the following. 

The map of Gran Sasso overburden as a function of direction and the description of its 



chemical composition are reported in Ref. l22fl . The detector simulation is based on the 



CERN package GEANT The folding of simulated events with the detector simulation 



is performed according to a variance reduction method |23| to minimize statistical losses and 
reducing possible systematic errors. 

We generated 3.6- 10 8 primary interactions in the total energy range 3-hlO 5 TeV, assuming 
the "MACRO-fit" primary mass composition model H|l(|, in which five mass groups (p, He, 
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CNO, Mg and Fe) are considered. Simulated data are produced with the same format as 
real data and are processed using the same analysis tools. After the application of the same 
cuts as for real data, a sample of about 7.0- 10 5 muon pairs survived, corresponding to about 
645 days of MACRO live time. 

In Fig. |] the comparison between the experimental and simulated decoherence curve 
inside the detector is shown. Curves are normalized to the peak of the dN/ dD distribution. 
The remarkable consistency of the two curves demonstrates the HEMAS code capability 
to reproduce the observed data up to a maximum distance of 70 m. The bump in the 
experimental distribution around 40 m is due to the detector acceptance and is visible also 
in the simulated data, thus confirming the accuracy of our detector simulation. We also 
notice that, despite the application of cut C4, there is a non- negligible discrepancy between 
the experimental and simulated data in the first two bins of the distribution of (34 ±2)% 
and (10±1)%, respectively. Such a discrepancy is not predicted by any model, since at short 
distances, apart from detector effects, the shape of decoherence distribution is dictated by the 
solid angle scaling: dN / dD 2 \u^Q ~ const., while the relevant properties of the interactions 
under investigation manifest themselves in the shape at large distances. The origin of this 
discrepancy will be discussed in detail in the last section of this paper, where other sources 
of contamination in the real data sample will be taken into account. 



IV. UNFOLDING PROCEDURE 



The agreement of the Monte Carlo and data shown in Fig. |4| proves that the simula- 
tion is consistent with observation and that the detector structure is well reproduced. A 
detector-independent analysis is required in order to subtract the geometric effects peculiar 
to MACRO, and allows a more direct comparison with other analyses and/or hadronic in- 
teraction models. This is accomplished by a correction method, built with the help of the 
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Monte Carlo simulation, to unfold the "true" decoherence function from the measured one 
in which geometrical containment and track reconstruction efficiencies are considered. 

In the previous decoherence studies the unfolding procedure was based on the 

evaluation of the detection efficiency for di-muon events generated by the Monte Carlo 
with a given angle and separation. Although this method is composition independent and 
allowed us to determine the detector acceptance with high statistical accuracy, it introduces 
systematic effects that have so far been neglected. In particular, in a multi-muon event it 
may happen that in a given projective view and in a particular geometrical configuration 
one muon track is "shadowed" by another. To avoid this effect, we adopt the following new 
unfolding method: the efficiency evaluation is performed considering the whole sample of 
events generated with their multiplicities. For a given bin of (D,9,4>), where D is the muon 
pair separation and (6,<p) is the arrival direction of the event, we calculate the ratio 



between the number of pairs surviving the selection cuts N tn and the number of pairs 
inserted in the detector simulator N out . In principle, this choice of e could be dependent 



dependent on the muon density and hence on its multiplicity, which in turn is correlated with 
the average atomic mass (A) of the primary. To check the systematic uncertainty related 
to this possibility, we evaluated the decoherence distributions obtained by unfolding the 
experimental data assuming the "heavy" and "light" composition models. Fig. [5] shows the 
relative comparison of the shape of the unfolding efficiencies as a function of pair separation, 
integrated in (cosO, 0) after the normalization to the peak value. In the same plot we present 
the unfolding efficiency calculated with the method used in Ref . |7||| . Considering the effect 
of the normalization, we observe that this method tends to overestimate the efficiencies in 
the low distance range, a consequence of the shadowing effect as explained in Section 2. 





on the primary mass composition model, since for a fixed distance D the efficiency 
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The unfolded decoherence is given by 

(dN\ N*"(D,e,<i>) f4 2) 

\ dD ),nf < D ^® ' 

where N exp (D, 9, (f>) is the number of muon pairs detected with a separation D. In practice, 
we used 50 windows in (cos9, (f>) space (5 and 10 equal intervals for cos 9 and <fi, respectively). 



The ability to evaluate the integral |Q| for separate and independent windows constitutes 
a powerful check of the systematics related to the decoherence dependence on the variables 
(cos9,(f>). Unfortunately this is not possible for r larger than 45 m, due to insufficient 
statistics. In that case the observables j\f m ^]\f out an d N exp are integrated over (cos9,4>). 
We verified that the systematic error introduced by that choice is smaller than the present 
statistical error in that distance range. 

Finally, unfolded experimental data obtained with the MACRO-fit model are directly 
compared with the Monte Carlo simulation (Fig. |6]). The two curves are in good agreement 
although the disparity in the first bin of the distribution remains unresolved (see Section 
VI). The experimental values of the dN/dD distributions, normalized to the peak value, are 
reported in Table |. 



V. UNCERTAINTIES OF THE HADRONIC INTERACTION MODEL 

The present work, as are others from MACRO, is extensively based on the HEMAS 
code. This was explicitly designed to provide a fast tool for production of high energy 
muons (E^ >500 GeV). However, as mentioned before, the interaction model of HEMAS is 
based on parametrizations of existing accelerator data and therefore is subject to the same 
risks of all this class of simulation codes. In particular, important correlations might be 
lost, or wrong, or the necessary extrapolations required by the specific kinematic regions of 
cosmic ray physics could yield unrealistic results. This remains a central problem of cosmic 
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ray physics. For this reason in the last few years general interest has grown in "physically 
inspired" simulations. These are based upon theoretical and phenomenological models like 
QCD and the Dual Parton Model capable of properly constraining the predictions where 
data do not exist, without the introduction of a large number of free parameters. It is worth- 



while to mention the attempt to merge the DPMJET model p5| into the shower simulation 
of HEMAS |[26|| , and the interface of the CORSIKA shower code with different mod- 
els, like HDPM (the original interaction model of CORSIKA), VENUS |28|, QGSJET @, 
SIBYLL |[30| , and the afore mentioned DPMJET. A review of general results obtained using 
CORSIKA with those models has been provided by the Karlsruhe group [^TJ. A common 
feature of all these models is the more or less direct reference to the Regge-Gribov theories 



32| for the soft contribution (low P±). It must be stressed that such a phenomenological 



framework, by its nature, provides only predictions for the longitudinal properties of the in- 
teraction. The transverse structure leading to the specific P± distribution is not constrained 
by the theory, except for the higher P± phenomena, where perturbative QCD can be used 
(this is of small relevance in the primary energy region addressed by the MACRO data). 
Once again, the model builders have to be guided mostly by experimental data, introduc- 
ing a-priori functional forms along with their additional required parameters. Some of the 
quoted models introduce proper recipes for the continuity between the soft and perturbative 
QCD regimes, and also specific nuclear phenomena like the Cronin effect mentioned above 
(see for instance In practice, the only possibility to evaluate a systematic uncertainty 

associated with the simulation model (at least those concerning the transverse structure of 
the showers) is to compare the predictions from all these models, HEMAS included. For this 
purpose, since the Karlsruhe report did not address this point, we have performed test 
runs with some of the models interfaced to CORSIKA, to which PROPMU p| has also been 
interfaced by us for muon transport in the rock overburden. A full simulation with all the 
other codes was outside our present capability, so we limited ourselves to comparisons at a 
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few fixed primary energies, and at fixed primary angles of 30° in zenith and 190° in azimuth. 
These correspond to an average rock overburden of ~ 3200 hg/cm 2 . In Tables [II] we show 
this comparison for a few representative average quantities for 3 different primary proton 
energies. We have considered the average depth of the first interaction X, (P±) for pions 
coming from the first interaction, the average production slant height of muons surviving 
underground (the decay height of their parent mesonsB.), the average distance of the muons 
from shower axis (R) and the average underground decoherence (D). Before discussing the 
results, it is important to remark that as far single interactions are concerned, all the models 
considered give a P± distribution following, with good approximation, the typical power law 
suggested by accelerator data. This is oc 1/(P± + Po) a , although with somewhat different 
parameters for different models. Older models, like those predicting a simple exponential 
distribution for P±, cannot reproduce the muon lateral distribution observed in MACRO 
data 0. 

In the energy range of 100-1000 TeV, to which most of MACRO data belong, the resulting 
differences in the average muon separation do not exceed 20%. These discrepancies seem 
to reduce at higher energy, while they appear much larger at few tens of TeV. DPMJET 
is probably the only model predicting a higher average separation than HEMAS. A precise 
analysis of the reasons leading to the differences among models is complicated. However, we 
note that HEMAS gives in general higher values of average P± than the other models. The 
only exception is indeed DPMJET, which, as mentioned before, pays particular attention 
to the reproduction of nuclear effects affecting the transverse momentum, as measured in 
heavy ion experiments |p3 . On the other hand, the effect of this large P± on the lateral 



distribution of muons is moderated in HEMAS by a deeper shower penetration (the inelastic 



*CORSIKA does not allow direct access to the production height of parent mesons, which would 
be more interesting for our purposes 
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cross section is based on Ref. ||34|| ); in general HEMAS exhibits a somewhat smaller height 
of meson production. 

Similar features in the comparison of models are also obtained for nuclear projectiles. It 
is therefore conceivable that, for the same primary spectrum and composition, not all the 
models considered could reproduce the MACRO decoherence curve. Thus the best fit for 
spectrum and composition as derived from the analysis of muon multiplicity distribution in 
MACRO will also probably differ according to the model. 

At least in part, the decoherence analysis can disentangle different ranges of longitudinal 
components of the interaction from the transverse ones, if this is performed in different zenith 
angle and rock depth windows. In fact, larger zenith angles correspond (on average) to larger 
muon production slant heights. This is a geometrical effect due to the greater distance from 
the primary interaction point to the detector for large zenith angle and consequently to the 
greater spreading of the muon bundle before reaching the apparatus. Larger rock depths 
select higher energy muons and consequently higher average energy of their parent mesons. 
The average separation decreases with the rock depth since, qualitatively, the longitudinal 
momentum (Pm) increases linearly with energy while (-Pi) increases only logarithmically. 
The overall result of increasing rock depth is the production of final states in a narrower 
forward cone, decreasing the muon pair average separation observed at the detector level. 

In Fig. [7] and |3| the unfolded decoherence function is compared to the HEMAS prediction 
for different zenith and rock depth intervals. In Table PH the average separation (D) is 
reported as a function of cosO and rock depth for fixed rock depth and zenith, respectively. 
In the same table we report the average values of slant height of first interaction (X), muon 
production slant height (Hfj), energy (E p ) and transverse momentum (P±) of the parent 
mesons, as obtained from the HEMAS Monte Carlo in the same windows. 

The agreement between the results and the Monte Carlo in separate variable intervals 
reinforces our confidence in the capability of HEMAS to reproduce the significant features 
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of shower development. This also allows us to exclude the existence of significant systematic 
errors related to this analysis. 



VI. THE CONTRIBUTION OF THE fj,± + N -> p± + N + /i + + /z~ PROCESS AT 

SMALL DISTANCES 

The capability of the MACRO detector to resolve very closely spaced tracks permits the 
extension of the decoherence analysis to a distance region hardly studied in the past. The 
mismatch between experimental and simulated data in this region (D < 160 cm) has been 
emphasized earlier in our discussion. In Section 2 a solution was attempted, permitting us 
to discard, with high efficiency, those tracks originating from secondary particle production. 
However, Fig. |4] and Fig. |6] show that other sources of contamination in the first bin of the 
decoherence function are responsible for the discrepancy. 

The process of muon pair production by muons in the rock, y u ± + N — > ^ + iV + fi + + \i~ , 



is a natural candidate. As pointed out in [16|, at the typical muon energy involved in 



underground analyses (E^ ~ 1 TeV) and for very large energy transfer, the cross section for 
this process is non-negligible with respect to e + e~ pair production. An analytic expression 
for the muon pair production cross section is given in [0,0 . In order to test the hypothesis, 
such a cross section has been included in the muon transport code PROPMU. Assuming a 
muon flux with energy spectrum E~ 3 ' 7 and minimum muon energy E™ m = 1.2 TeV at the 
surface, and considering the actual mountain profile, we generated a sample of 10 7 muons 
corresponding to 3666 h of live time. About ~ 3.0 • 10 6 muons survived to the MACRO level, 
5360 of which were generated by muon pair production processes. The average separation 
of these muon pairs is (128 ± 1) cm, and their average residual energies are (657 ± 14) GeV 
and (145 ± 3) GeV, respectively, for the main muon and the secondary muon samples. We 
propagated the muons surviving to the MACRO level through the GEANT simulation and 
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we applied the same cuts specified in Section 2. Finally, the number of events was normalized 
to the live time of real data. 

In Table |V| we report the number of weighted muon pairs in the first bins of the experi- 
mental and simulated decoherence distributions (in the form dN/dD). The effect of standard 
cuts, of cut C4, and of the subtraction of the muon pair production process are shown in 
order. In each case, we indicated in percentage the bin populations with respect to the peak 
of the distribution and the discrepancy with respect to the Monte Carlo predictions. 

In Fig. [5] we compare the simulated decoherence curve with the data corrected for the 
muon pair production effect. Despite the approximation introduced in our test, it seems that 
the proposed muon pair production process can account for most of the observed discrepancy 
in the low distance range. This is also shown in the inset of Fig. ||| where the distribution 
of relative distance for the muon pairs in excess of the data (after subtraction of HEMAS 
prediction) is compared to the expectation from simulated muon pair production. 

An excess at small pair separation is also predicted in exotic processes, like multi-W 



production by AGN z/s, as suggested in |36|]. However, according to this reference, muons 
from W— decay have an average energy of ~80 TeV. These muons would survive 
underground with a residual energy much higher than that of standard muons, producing 
local catastrophic interaction in the detector, making difficult their identification as a pair. 
On the contrary, the explanation proposed here is based on a pure QED process that does 
not require any additional physics. 



VII. DISCUSSION AND CONCLUSIONS 



We have obtained an improved experimental underground decoherence function using 
high energy muons (E^ > 1.3 TeV) up to a maximum distance of about 70 m. It is hard to 
conceive in the near future a large area underground experiment capable of improving the 
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sensitivity reached in this decoherence study. 

A new unfolding of the experimental distribution confirms the results obtained with the 
analysis within the detector. 

The ability to resolve closely spaced muon tracks allows an investigation of the deco- 
herence function at small separations. Apart from the negligible contamination of hadro- 
production by muons (which will be the subject of a future work), we found that a relevant 
contribution is made by the process /i ± + iV — > /i ± + iV + yU + + //~. The inclusion of this 
interaction in the simulation reproduces, in both a qualitative and a quantitative way, the 
experimental data. 

The agreement of the overall distribution shape for experimental and simulated data from 
HEMAS is excellent. The possible excess at high muon separations suggested in the previous, 
preliminary, analyses ]7||| was due to an imperfect unfolding procedure, and is now excluded. 
These results both in the integrated distribution and in those from separate intervals of zenith 
and rock depth, shows that HEMAS gives a reasonable account of the cascade development 
and that it is not necessary to introduce any anomalous P± production in the Monte Carlo 
to reproduce these data. However, the other interaction models considered for comparison, 
while reproducing similar behaviour, in general give different combinations of transverse 
momentum and production height. Discrimination among the different models may be 
possible only after a complete simulation and analysis of MACRO data with each of the 
codes. Therefore the present work, representing a final data reduction and analysis, provides 
a valuable benchmark for future analysis dedicated to the investigation of the properties of 
high energy interactions and to the evaluation of different shower models in the primary 
energy region spanning from a few tens to a few thousands TeV/nucleon. The detector 
independent analysis described here will make this task easier. 



20 



REFERENCES 



[1] T.K. Gaisser, Proc. Vulcano Workshop 1992 "Frontier Objects in Astrophys. and Particle 
Phys.", 40, 433 (1993). 

[2] P. Sokolsky, Introduction to Ultra High Energy Cosmic Ray Physics, Addison & Wesley 
Publishing Co. (1988), p. 131. 

[3] C. Forti et al, Phys. Rev. D42 3668 (1990). 

[4] H.E. Bergeson et al., Phys. Rev. Lett. 35 1681 (1975). 

[5] G. Battistoni et al, LNGS-95-09 (1995); G. Battistoni et al, in Proceedings of the XXIV 
Int. Cosmic Ray Conf, Roma, Italy, 1995, ed. N. Iucci et al., Arti Grafiche Editoriali, 
Urbino, (1995), Vol. 1, p. 508. 

[6] MACRO Collaboration, S.P Ahlen et al., Nucl. Instr. Meth. A324 337 (1993). 

[7] MACRO Collaboration, S.P. Ahlen et al., Phys. Rev. D46 4836 (1992). 

[8] MACRO Collaboration, S.P. Ahlen et al., Proceedings of the XXIII Int. Cosmic Ray Conf., 
Calgary, Canada, 1993, ed. R.B. Hicks et al., World Scientific, Singapore, (1994), Vol. 2, 
p. 93. 

[9] MACRO Collaboration, M. Ambrosio et al., Phys. Rev. D56 1407 (1997). 

[10] MACRO Collaboration, M. Ambrosio et al., Phys. Rev. D56 1418 (1997). 

[11] G. Auriemma et al., Proceedings of the XXI Int. Cosmic Ray Conf., Adelaide, Australia, 
1989, ed. R.J. Protheroe, Graphic Services, Northfield, Australia, (1990), Vol. 9, p. 362. 

[12] J.W. Cronin et al., Phys. Rev. Dll 3105 (1975). 

[13] M. Bardadin-Otwinowoka et al., Proceedings of the Vienna Conference on Elementary 
Particles, Vol. 1, 628 (1963). 

21 



[14] S0UDAN2 Collaboration, S.M.S. Kasahara et al., Phys. Rev. D55 5282 (1997). 

[15] MACRO Collaboration, M. Ambrosio et al., Proceedings of the XXV Int. Cosmic Ray 
Conf., Durban, South Africa, 1997, ed. M.S. Potgieter et al., Wesprint, Potchefstroom, 
(1997), Vol. 6, p. 357. 

[16] V.A. Kudryavtsev and O.G. Ryazhskaya, INFN-AE-97-08 (1997); V.A. Kudryavtsev and 
O.G. Ryazhskaya, in Proceedings of the XXV Int. Cosmic Ray Conf., (1997), Vol. 6, p. 
405. 

[17] A. Fasso et al., An update about FLUKA, Proc. 2nd workshop on Simulating Accelerator 
Radiation Environment, SARE-2, CERN-Geneva, 9-11, October, 1995. 

[18] R. Brun et al., CERN GEANT3 User's guide, DD/EE/84-1 (1992). 

[19] P. Lipari and T. Stanev, Phys. Rev D44 3543 (1991). 

[20] J. Engels et al., Phys. Rev. D 46, 5013 (1992). 

[21] G. Battistoni et al., Nucl. Instr. Meth. A394 136 (1997). 

[22] MACRO Collaboration, M. Ambrosio et al., Phys. Rev. D52 3793 (1995). 

[23] G. Battistoni et al., Astroparticle Phys. 7 101 (1997). 

[24] A. Capella et al., Phys. Rep. 236 225 (1994). 

[25] J. Ranft, Phys. Rev. D51 64 (1995). 

[26] G. Battistoni et al., Astropart. Phys. 3 157 (1995). 

[27] J.N. Capdevielle et al., the Karlsruhe extensive air shower simulation code CORSIKA, 
KFK Report 4998 (1992). 

[28] K. Werner, Phys. Rep. 232 87 (1993). 



22 



[29] N.N. Kalmikov et al, Physics of Atomic Nuclei 58 1728 (1995). 
[30] R.S. Fletcher et al, Phys. Rev. D50 5710 (1994). 

[31] J. Knapp, D. Heck and G. Schatz, Comparison of hadronic interaction models used in Air- 
shower simulations and their influence on shower development and observables, Karlsruhe 
report FZKA 5828 (1996). 

[32] V.N. Gribov, Sov. Phys. JETP 26 414 (1968); V.N. Gribov and A. A. Migdal, Sov. J. 
Nucl. Phys. 8 583 (1969); see also M. Baker and K.A. Ter-Martirosyan Gribov's Reggeon 
calculus: its physical basis and implication, Phys. Rep. 28C 1 (1976). 

[33] See the review of H.R Schmidt and J. Schukraft, GSI-preprint GSI-92-19 (1992). 

[34] L. Durand and H. Pi, Phys. Rev. Lett. 58 303 (1987). 

[35] S.R. Kelner, Sov. J. Nucl. Phys. 5 778 (1967). 

[36] D.A.Morris and A.Ringwald, Astropart. Phys. 2 43 (1994). 



23 



FIGURES 




E 

51 60 



50 
40 




(b) 



30 



20 



10 



40 50 
Distance (m) 

FIG. 1. Average separation of underground muon pairs at Gran Sasso depth, as a function 
of (P±) of the parent mesons (a) and of the slant height in the atmosphere (b). The results are 
obtained with the HEMAS Monte Carlo. 
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FIG. 2. The change in the experimental decoherence function induced by cut C4- The data 
indicate the fractional deviation between the experimental decoherence function before and after the 
application of the cut. 
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FIG. 3. Percentage of reconstructed real events (white points) and unambiguously associated 
muon pairs (black points) as a function of event multiplicity. 
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FIG. 4. Experimental (black points) and simulated (white points) decoherence function, nor- 
malized to the peak of the dN/dD distribution. The second to last points of the two distributions 
coincide. 
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FIG. 5. Comparison of the unfolding efficiencies as a function of pair separation for different 
composition models. The curves are normalized to the peak value. For comparison, we include the 
efficiency evaluated with the method used in previous analyses (white points). 
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FIG. 6. Unfolded experimental decoherence distribution compared with the infinite- detector 
Monte Carlo expectation, computed with the HEMAS interaction code and the MACRO-fit pri- 
mary composition model. Black squares represent data above 45 m (integral form unfolding). 
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FIG. 7. Unfolded decoherence functions compared with Monte Carlo simulations for different 
cosO windows. The vertical scale is in arbitrary units. 
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FIG. 8. Unfolded decoherence functions compared with Monte Carlo simulations for different 
rock depth windows. 
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FIG. 9. The low distance region of the experimental decoherence function, before and after 
the subtraction of the secondary muon sample, and comparison with the Monte Carlo simulation. 
The inset shows the distribution of relative distance for muon pairs in excess of the data after the 
subtraction of the HEMAS prediction, as compared to the expectation from simulated muon pair 
production in the rock. 
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TABLES 



D (cm) Error 

80 0.4254 0.1034 10" 1 

240 0.8855 0.1591 10" 1 

400 1.0000 0.1738 10" 1 

560 0.9591 0.1719 10" 1 

720 0.8148 0.1554 10" 1 

880 0.6730 0.1419 10" 1 

1040 0.5595 0.1350 10" 1 

1200 0.4341 0.1216 10" 1 

1360 0.3410 0.1118 10" 1 

1520 0.2939 0.1144 10" 1 

1680 0.2198 0.4950 10~ 2 

1840 0.1828 0.4626 10~ 2 

2000 0.1578 0.4476 10~ 2 

2160 0.1283 0.4134 10~ 2 

2320 0.1047 0.3906 10~ 2 

2480 0.9348 lO" 1 0.3881 10~ 2 

2640 0.7438 10" 1 0.3599 10~ 2 

2800 0.5847 10" 1 0.3211 10~ 2 

2960 0.5168 10" 1 0.3238 10~ 2 

3120 0.4173 lO" 1 0.2994 10~ 2 

3280 0.4113 10" 1 0.3125 10~ 2 

3440 0.2582 lO" 1 0.2585 10~ 2 

3600 0.2315 lO" 1 0.2444 10~ 2 

3760 0.2260 10" 1 0.2522 10~ 2 
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3920 


0.2004 lO" 1 


0.2501 10~ 2 


4080 


0.1289 10- 1 


0.2023 10~ 2 


4240 


0.1419 lO" 1 


0.2319 10~ 2 


4400 


0.9105 10~ 2 


0.1775 10~ 2 


4560 


0.6776 10~ 2 


0.1597 10~ 2 


4720 


0.3080 10~ 2 


0.1035 10~ 2 


4880 


0.2990 10~ 2 


0.1121 10~ 2 


4640 


0.1128 10- 1 


0.3340 10~ 2 


4960 


0.8697 10~ 2 


0.4549 10~ 2 


5280 


0.7108 10~ 2 


0.4098 10~ 2 


5600 


0.3689 10~ 2 


0.1534 10~ 2 


5920 


0.4190 10" 2 


0.2742 10~ 2 


6240 


0.2991 10~ 2 


0.2703 10~ 2 


6560 


0.7539 10~ 3 


0.9004 10-3 


TABLE I. 


Tabulation of the the unfolded decoherence 


distribution as measured by MACRO. 



The data points are normalized to the point of maximum. 



p-Air, 20 TeV 



Code 


(X first) 


(P±) ^ 




(R) 


(D) 




(g/cm 2 ) 


(GeV/c) 


(km) 


(m) 


(m) 


HEMAS 


51.4 


0.40 


24.1 


7.9 


12.7 


CORSIKA/DPMJET 


44.4 


0.42 


25.6 


10.1 


13.9 


CORSIKA/QGSJET 


45.7 


0.39 


24.3 


7.3 


10.0 


CORSIKA/ VENUS 


48.3 


0.35 


24.5 


7.4 


8.3 


CORSIKA/SIBYLL 


50.9 


0.37 


23.5 


7.2 


11.5 


p-Air, 200 TeV 


Code 


(A first) 


(P±) *± 


w 


(R) 


(D) 




(g/cm 2 ) 


(GeV/c) 


(km) 


(m) 


(m) 


HEMAS 


56.1 


0.44 


20.6 


5.3 


8.0 


CORSIKA/DPMJET 


53.9 


0.43 


21.7 


6.2 


8.8 


CORSIKA/QGSJET 


52.8 


0.41 


21.4 


5.5 


7.8 


CORSIKA/VENUS 


60.2 


0.36 


20.9 


5.3 


7.5 


CORSIKA/SIBYLL 


55.2 


0.41 


20.2 


5.2 


7.3 


p-Air, 2000 TeV 


Code 


(A fi rs t) 


(P±) *± 




(R) 


(D) 




(g/cm 2 ) 


(GeV/c) 


(km) 


(m) 


(m) 


HEMAS 


63.0 


0.50 


16.3 


4.1 


6.0 


CORSIKA/DPMJET 


60.0 


0.42 


18.5 


4.9 


6.4 


CORSIKA/QGSJET 


63.1 


0.44 


17.7 


4.2 


5.6 


CORSIKA/VENUS 


66.7 


0.36 


16.8 


4.1 


5.3 
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CORSIKA/SIBYLL 60.3 0.44 17.0 4.4 5.6 

TABLE II. Comparison of a few relevant quantities concerning the lateral distribution of un- 
derground muons at the depth of 3200 hg/cm 2 , from proton primaries at 20, 200 and 2000 TeV, 
30° zenith angle. The statistical errors are smaller than the last reported digit. 
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3750<h<4150 (hg/cm 2 ) 







.5<cos6><.6 


.6<cos0<.7 


.7<cos6><.8 


.8<cos6><.9 


.9<cos0<l. 


EXP 


(D) (m) 


13.2 ±2.3 


11.4 ±2.2 


10.3 ±2.2 


8.5 ± 1.9 


7.5 ±1.9 




(D) (m) 


12.8 ± 1.4 


12.0 ± 1.3 


10.1 ± 1.2 


8.8 ±1.2 


7.8 ±1.1 




(X) (km) 


65.9 ±0.2 


57.0 ±0.2 


51.2 ±0.3 


42.5 ±0.4 


37.0 ±0.5 


MC 


W (km) 


41.6 ±0.3 


34.3 ± 0.3 


28.1 ±0.3 


23.8 ±0.3 


20.4 ±0.3 




(E p ) (TeV) 


4.1 ±0.2 


4.0 ±0.2 


4.0 ± 02 


3.9 ±0.1 


3.9 ±0.2 




(Pl) (GeV/c) 


0.56 ±0.01 


0.59 ±0.01 


0.57 ±0.01 


0.57 ±0.02 


0.57 ±0.01 



(a) 



0.8<cos#<0.9 







3350<h<3750 


3750<h<4150 


4150<h<4550 


4550<h<4950 






(hg/cm 2 ) 


(hg/cm 2 ) 


(hg/cm 2 ) 


(hg/cm 2 ) 


EXP 


(D) (m) 


9.4 ±2.1 


8.5 ± 1.9 


7.3 ± 1.6 


6.2 ± 1.6 




(D) (m) 


9.7 ±3.4 


8.8 ± 1.2 


7.7 ±1.1 


7.1 ± 1.1 




(X) (km) 


42.7 ±0.4 


42.5 ± .4 


45.9 ±0.3 


43.76 ± 0.3 


MC 


W (km) 


23.7 ±0.3 


23.8 ±0.3 


24.6 ± 0.3 


25.1 ±0.5 




(E p ) (TeV) 


3.6 ±0.1 


3.9 ±0.1 


4.4 ±0.1 


4.8 ±0.02 




(Pl) (GeV/c) 


0.56 ±0.02 


0.57 ±0.02 


0.58 ±0.02 


0.58 ± 0.02 



(b) 

TABLE III. Average separation between muon pairs (D) (in m) as a function of cosO (a) 
and rock depth (b). In each table the experimental data are compared to the expectations from the 
HEM AS Monte Carlo. For the same simulations, other averages of relevant quantities are reported. 
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0-80 cm 


80-160 cm 


160-240 cm 


240-320 cm 


320-400 cm 


PVyn Data 

J_J j\. JJ . 1 J CI L CI 




12491 


1 7569 


20514 




MP Data 


51 54 


9141 7 

Zj 1t:1 / 


33573 


40367 

1UOU 1 


49679 


Discrepancy after 


(55±2)% 


(16±2)% 

\ J- W 1 id / / \J 


(6±1)% 


(4±1)% 




noi iiiaiizaiioii 












F,xr> Data 4- CA 


3612 


11128 


1 6535 


1 9597 


19977 


MC Data + CA 


4848 


20346 


31932 


38425 


40660 




^o^znz j /o 


^iuin^ ^ /o 


yJZilZj J /o 






normalization 












Exp. Data + C4 + 


2193 


9264 


15462 


19190 


19842 


/i pair subtraction 












MC Data + C4 


4848 


20346 


31932 


38425 


40660 


Discrepancy after 


(8±7)% 


(7±3)% 


(0±2)% 


(2±2)% 




normalization 













TABLE IV. Number of weighted muon pairs in the first few bins of the experimental and simu- 
lated decoherence distributions. The discrepancy is the percentage difference between experimental 
and Monte Carlo values, normalized to the distribution maximum (last column). 
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